%Generates instances of autoregressive processes with 
%sparsity patterns in the inverse spectrum which correspond to 
%the two step graph generated form the graph defined by the matrix G.
%Outputs an N time ordered samples of a realization in a matrix x of N by n,
%outputs the auto regressive coefficients A = [I -A_1,..,-A_p], O=Omega
%Equivalent model for white noise innovations B = [B0 B_1,...,B_p] where B0*B0 = Omega
%The AR order p, the variates dimension n.
%The windowed estimate of the covariance matrix C and saves it to a file at 'file'

%clear all
%Set the random seed
%rand_stream = RandStream('mt19937ar','Seed',10);
%RandStream.setDefaultStream(rand_stream);
%Generate a graph by forming a cycle of n vertices, 
%introduce self edges and edges to the nodes one step ahead

if ~exist('G') ||  ~exist('N') || ~exist('n') || ~exist('p')
  	N = 500;
	n = 6;
	p = 4;
  	fprintf('Will use a cycle of length %i as sparsity pattern\n',n);
	G = eye(n) + diag(ones(n-1,1),1) +diag(ones(n-1,1),1);
	G(1,n) = 1;
	G(n,1) = 1;
	G = sparse(G);
else
  fprintf('Will use the sparsity pattern of the variable G in scope\n')
  if size(G,1)~=size(G,2), fprintf('G is not square\n'); break; end
  if ~issparse(G), fprintf('G is not sparse\n'); break; end
  if n ~= size(G,1), fprintf('n does not correspond with the size of G will change n  to %i\n',size(G,1));n = size(G,1) ;end
end
  Gsq = G'*G;
if nnz(Gsq) == n*n, fprintf('Warning G^2 is not sparse!\n'); end

stable = false;
%Iterate untill we get a stable model
while ~stable

	%Generate p+1 random sparse matrices with
	%sparsity pattern equal to that of G and
	%store them into a matrix formed by the
	%concatenation of all the matrices
	B = spalloc(n,n*(p+1),n*n*(p+1));
	
	%Generate B0, this matrix has to be full rank
	%we will add a small multiple of the identity,
	%this will affect the sparsity of the end spectrum.
	Vs = 0.1*sprandsym(G+G'); 
    B(:,1:n) =  Vs + diag(0.5+sum(abs(Vs)));
    
	%Generate the rest of the parameters
	for i = 1:p
	 B(:,(i)*n+1:(i)*n+n) = 1/p*sprandn(G);
	end		
	B0 = B(:,1:n);
	iB0 = inv(full(B0));

	%Generate the regression coefficients
	A  = iB0*B(:,n+1:end);

	%Form the companion matrix 
	Comp = [-A;[eye(n*(p-1)),zeros(n*(p-1),n)]];
	ei= eigs(Comp); 
	me = max(abs(ei)); 
	fprintf('Maximum eigenvalue of the companion matrix %d\n',me)
	if max(abs(ei)) < 1-eps
		stable = true;
		fprintf('Selected stable model with maximum eigenvalue %d\n',me);
	end
end


%Generate the Ar process
Xt = zeros(p*n,1);


burn = ceil(log(eps)/log(me));
fprintf('Burn in for %d iterates\n',burn) 
%Estimate the length of the burn in 
%Burn in process
X = [];
for j = 1:burn
	Xn = -A*Xt+iB0*randn(n,1);
	Xt = circshift(Xt,n);
	Xt(1:n) = Xn;
end
fprintf('Generating realization \n')
%Generate the realization with the B coefficients
for j   = 1:N
	X   =  [X,Xt((p-1)*n+1:end)];
	Xn  = -A*Xt + iB0*randn(n,1);
	Xt  = circshift(Xt,n);
	Xt(1:n) = Xn;	
end


%Form the coefficient matrices
A = [eye(n) A];
%B = [B0 B];
%free the symbol X
x = X';

%compute the sparsity pattern 
Y = compute_D(B'*B,n,p);
Sp = abs(Y(:,1:n));
for i = 1:p, Sp = Sp + abs(Y(:,n*i+1:n*i+n)); end;
Sp = Sp ~= 0;


%Generate the windowed estimate of the covariance matrix
[N,n] = size(x);
C = windowed_est(x, p);
%make sure numerical error does not make it non symmetric
C = (C+C')/2;
%Build the sparsity pattern

clear Comp;
clear W;
clear burn;
clear X;
clear Xt;
clear Xn;
clear B0;
clear H;
clear stable;
clear me;
clear i;
clear j;

file = '../Data/graph_7_4_50.mat';
save(file)
fprintf('Results saved to %s\n',file)
%no need to do the p = k change, data is now lower case x with column per variable and row per sample.


